Plastic flow in a sheared polycrystalline solid using phase field model 
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Plastic deformation in solids induced by external shear stress is of huge practical interest. Presence of local 
crystalline order in polycrystals, consisting of many grains, distinguishes its deformation pattern from that of 
amorphous materials. Despite strong anisotropy, induced by external stress, the plastic flow and the consequent 
deformation field show strong dynamical heterogeneity. The distribution P(u) of particle displacements (u) 
shows three distinct regimes including a power law scaling regime at moderate displacements. Using a phase 
field simulation we show how polycrystals generate saddle and vortex like flow patterns, which hitherto have 
been termed as elementary plastic events in the context of amorphous materials. Interestingly, such events 
here find natural explanation in terms of the underlying dislocation dynamics. We also characterize the spatial 
distribution of the flow field using Okubo- Weiss measure. 
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Solids undergo plastic deformation and subsequently flow 
when they are subjected to stresses beyond their elastic limit. 
Understanding plastic flow in polycrystalline solids has huge 
practical importance, starting from applications in material 
processing to understanding earthquake dynamics. Metals are 
deformed routinely in various industrial applications, for ex- 
ample, by rolling or pressing them into sheets. In earthquake 
dynamics, it is the shear stress generated by the movement of 
tectonic plates which causes plastic flow. Such deformations 
and flows under constant shear stress as well as constant shear 
strain rates have been widely studied both by lab experiments 
HI] and simulations J2l- Experiments on crystalline ice has es- 
tablished that acoustic emission which is a signature of dislo- 
cation avalanches in stressed solid is strongly intermittent and 
shows power law distribution of the energy released by such 
activity. Simulations using dislocation dynamics model and 
crystal phase field models have reproduced this JH 01. Simu- 
lation of amorphous solids ®^ has also shown signature of 
spatially intermittent behavior which often organize into large 
flows in the form of vortices and saddles. 

Plastic flow in a polycrystal ranges over multiple length 
scales starting from the motion (glide/climb) of single disloca- 
tions at the nanometer scale to movement of grain boundaries 
at micron and larger scales. Also they are highly anisotropic 
biased by the direction of the macroscopically applied exter- 
nal stresses. We are interested in the nonequilibrium dynamics 
of the plastic flow which is generated when the solid is sub- 
jected to external shear stress. While most studies have dis- 
cussed the displacement field which results from a quasi-static 
strain we subject the solid to a constant strain rate (01. 

We employ phase field crystal (PFC) model to simulate a 
sheared polycrystalline solid in two dimensions. The strength 
of the phenomenological PFC model is that it can study 
dynamics of crystalline solids at the microscopic (atomic) 
length scales but diffusive time scales (much longer than vi- 
brational time scales). Also here dislocations are generated 
spontaneously without any ad-hoc rules being imposed. PFC 



model have been augmented (modified PFC) with an accel- 
eration term in the damped Langevin equation. This reintro- 
duces the fast acoustic waves into the dynamics. PFC and 
modified-PFC have been successful in reproducing impor- 
tant phenomenology of grain-boundary energy dUl, pre melt- 
ing transition jjj], dislocation glide iflOll to mention a few. It 
has also been applied to stud y g lass transition time scales 111 ill 
and liquid crystal dynamics 1 1 211 . Modified PFC has also been 
derived ]8[] from microscopic density functional theory. 

We confine the square shaped solid between two parallel 
plates, at y = and y = 2H which are moved at uniform 
speeds v^x and — vqx, respectively. Along the x direction we 
apply periodic boundary condition. To simulate the dynamics 
we use the modified phase field crystal (MPFC) model devel- 
oped in Ref [ 1 1,17] and follow a shearing scheme proposed by 
Chan et al. |2|]. The dimensionless equation of motion of the 
system is as follows : 



^ = a ^ K v 2 ^ + v 3 ) + Ky)^ + C ! (i) 



where ip is the conserved order parameter field, w(V 2 ) = 
r + (1 + V 2 ) 2 and ( is the conserved noise, a, (3 variables 
control the time scale of the phonon modes propagating in the 
solid and the degree of their damping. In this shearing scheme 
instead of moving only the top and the bottom surfaces of the 
solid, a drift velocity profile, v(y)x = vq exp(— y/X) for < 
y < H and v(y)x = t>o cxp(— (2H — 2/)/A) for H < y < 
2H, that decays exponentially away from the surface (towards 
the bulk) is applied on the solid. We hope that as long as 
the length scale of decay A <C H, the width of the solid, the 
physics in the bulk will not be affected by A. We examined the 
effect of this imposed drift by measuring the average velocity 
along the flow direction, \{v x (y))\ as a function of y (see 
inset of FigQ}c) and inferred that there exists a sizable bulk 
portion. All our studies focus on this bulk region. 

The results presented here are from simulations on a square 
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grid of size 256 x 256. No qualitative difference was found 
for a bigger grid (1024 x 1024), except that the data is much 
better averaged for the presented case due to reasonable run 
time. The parameters used are H x = H y = 256, a = 1, 
(3 = 1, A = 10, Vo = 0.3, r = -0.5, dx = 3tt/8, v = 0.45 
and dt = 0.025. Note that our X/H y is half of that in Ref|Q]. 
We used a pseudo-spectral scheme, combined with integrat- 
ing factor method, to solve the PDE in Eq. Q] We used two 
different initial conditions: a) a single crystal and b) a random 
distribution of grains of different sizes, grown by implanting 
many artificial nuclei in the supercooled liquid state. Under 
constant strain rate both led to the same (statistically) nonequi- 
librium, steady polycrystalline state with relatively smaller 
grains near the boundaries than in the interior (see Fig.l-a). 

In PFC model, the particles are identified as the minima 
of the scalar field xp(r, t). The grains have triangular lattice 
structure (see Fig[T|l-a), typical of 2D, with the local orienta- 
tion field given by the angle (9(r) G [0, 7r/3]. The maximum 
misorientation between two neighboring grains can be ir/Q. 
Local crystal orientation at each particle was obtained from 
the positions of its neighbors and keeping in mind the triangu- 
lar symmetry of the crystal. At grain boundaries crystal orien- 
tation changes. Using Delaunay triangulation, the number of 
nearest neighbors of each atom and the specific atoms linked 
to it was obtained. Dislocations were located by finding pairs 
of atoms with 5 and 7 neighbors on the Delaunay network. 

As mentioned in the introduction, collective motion in re- 
sponse to small shear strain applied at the boundaries of a solid 
I2H31, has been studied for amorphous material. But since the 
applied strain there JH-jil was quasistatic in nature (keeping 
the system always at equilibrium or at worst metastable equi- 
librium), the motion was athermal (i.e., T = 0). In contrast 
here we study the actual temporal dynamics in response to a 
constant strain rate. Collective motion of a many body sys- 
tem is often interpreted as the motion Under a constant strain 
rate, as ours, the system has no time to settle down in the lo- 
cal minima consistent with the global strain. The landscape 
changes at a much faster time scale compared to the relax- 
ational kinetics of the system. Thus the collective motion is 
nonequilibrium in nature and we focus on the steady state fea- 
tures. 

Several new features arise because of the polycrystalline 
nature of the sheared solid. The grains do not transmit hydro- 
dynamic stresses from one part of the system to the other like a 
fluid or even like an amorphous solids. A grain resists motion 
till the accumulated strain crosses its elastic limit and then 
either rotate with respect to its neighboring grains or break 
up into smaller grains. Therefore the scale of motion is con- 
trolled by the typical grain size, which becomes smaller at 
higher strain rates. Fig. [T}b shows the detailed velocity map 
of the particles in the bulk. Despite the strong bias along x 
(the shear direction) the flow field shows significant motion 
along y, giving rise to characteristic vortical and saddle type 
of motion as shown in Fig.[T]-b. Origin of such kind of motion 
will be discussed later. The inset of Fig. [TJ-c shows (v x (y)) 
as a function of y, which shows an emergent velocity profile 



in the bulk differing from the imposed drift velocity. Despite 
the shear induced anisotropy the flow is highly heterogeneous 
even within the bulk: very slow in the interior of the grains and 
fast at the grain boundaries. Since flows are driven by local 
stresses, high shear stresses are expected at the grain bound- 
aries which we will discuss next. The plastic flow also results 
in interesting displacement patterns which will be discussed 
later. 

Large scale molecular dynamics (MD) simulations on 
sheared amorphous material J^HH have demonstrated that 
plastic displacement patterns organize into large scale vortices 
when the externally imposed, global, shear strain is changed 
quasi-statically. Ref. |4j] has shown that an elementary, irre- 
versible, plastic event is a quadrupolar displacement pattern 
which has long range elastic effect. Picard et al 11411 . stud- 
ied the effect of a localized plastic event in a sheared 2D vis- 
coelastic medium. They analytically showed how, through 
elastic interaction, it induces a long range strain field which 
turns out to of quadrupolar nature. The event was assumed 
to be a localized plastic strain eijS(r) corresponding to a pure 
shear deformation, i.e., only the non-diagonal elements of 
were non-zero and equal. Generation of such a plastic strain 
matrix is at best plausible given the material is being globally 
sheared through its external boundaries. The MD simulations 
of Maloney and Lemaitre J51 01 on 2D sheared amorphous 
materials could nail down such local plastic events to charac- 
teristic motion of particles near a saddle with one stable and 
another unstable axis, and particles move towards and away, 
respectively, from the saddle along these axes. But it is still 
unclear how such singular points are created in the interior 
of the system due to global shear stress applied at its distant 
boundaries. 

We identify a microscopic mechanism which can dynam- 
ically give rise to quadrupolar patterns, essentially a saddle, 
in the displacement field of the particles. We show that such 
a pattern results from sideways approach of two oppositely 
charged edge dislocations towards each other. The sequence 
in Fig.|2]clearly shows time development of the displacement 
field leading towards a saddle as the dislocations approach 
each other. The saddle fades away after the dislocations an- 
nihilate (not shown here). Such long range quadrupolar pat- 
terns have been reported for noncrystalline (amorphous) ma- 
terial J3l also, but there being no dislocations in noncrystalline 
material microscopic origin of such a quadrupolar pattern re- 
mained unexplained. Quantitatively we know that the dis- 
placement field of an edge dislocation generates a displace- 
ment dipole where the positive and negative lobes are oriented 
along the axis connecting the atoms with coordination num- 
bers 5 and 7. Such a 5 — 7 pair is like a charge dipole and 
sideways arrangement of two such pairs form a quadrupole. 
Essentially these two dislocations have opposite burgers vec- 
tors b and —b along one of the symmetry axes of the crystal. 
In comparison, dislocations with the same burgers vector can 
line up in an "..5-7-5-7-5-7.." arrangement to form a disloca- 
tion wall (see Fig. QJa) which is rather stable. These walls 
are equivalent to high angle grain boundaries (Fig. QJa). The 
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FIG. 1: (Color online) (a) Spatial map of the local crystal orientation 8(r) with a range [0, tt/3]. Dislocations are indicated by light colors. 
They decorate the grain boundaries separating grains with high misorientation. Across grains with very low misorientation (cells) they are 
absent. This is consistent with the Frank condition (n oc sin 8) fl3h . which relates the line density of dislocations (n) along a grain boundary 
and the corresponding misorientation angle (9). (b) Velocity pattern of the particles in the bulk. Velocities near the moving boundaries are 
too large in magnitude to be represented in the same plot. Vortices are visible near the top and bottom boundaries while saddles (circled) are 
visible in the bulk. We omitted the velocity pattern near the boundary region since the velocity scale is much higher there compared to the 
bulk, (c) Shows shows semi-log plot of (v x (y)) as a function of y. To bring out the anti-symmetry of the average flow with respect to the 
mid-plane (y — H), we reflected (v x (y)) for y < H across the mid-plane and then used modulus (increasing y is directed towards the bulk. 
The solid line shows the imposed drift velocity v(y) ~ vo exp( — X/y). The upper inset shows the same data in a log-log plot which brings 
out the scaling behavior y~ a in the bulk. The lower inset shows the same data in a regular x — y plot. 




FIG. 2: (a) and (b) focus on a small region of the lattice. The circles represent the atoms and the dark arrows on the circles represent their 
velocity vectors while the light arrows represent the interpolated velocity field. The filled circles are the atoms with five or seven neighbors, 
indicating an edge dislocation, (a) shows initiation of a saddle as two oppositely charged dislocations approach, (b) Time development of the 
saddle as the dislocations get close to form a quadrupole. Note that 'b' has higher resolution than 'c'. (c) The superposed displacement field 
of two edge dislocations (with opposite Burgers vector) calculated using Eq. (|2j- The parameters used are b — 1 and v = 0.1. The dislocation 
positions are (0, 0) and (— 1, 0). 



quadrupolar structure discussed above can be quantitatively 
established by superimposing the elastic displacement fields 
of two dislocations located close by. Fig.[2]shows the resultant 
field from two dislocations located at f = (0, 0) and (—1, 0), 
where displacement field u(r) due to a dislocation at the ori- 
gin is given by lfl5ll 
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where v is the Poisson ratio and tan -1 f G [0, 2w]. It is worth 
mentioning that we have found saddle structures in the dis- 
placement profile even in places where there are no disloca- 
tions. Thus the above mechanism cannot be the only reason 
for saddles. 

In order to study the spatial distribution of vorticity and sad- 
dles (centres) in our 2D plastic flow we employ a technique 
borrowed from fluid turbulence [ iH 17|. For 2D inviscid, in- 
compressible flows the Okubo Weiss parameter is defined as 
A = dct(diVj). This is an invariant of the flow and can be 
recast as A = oj 2 — e 2 , where uj = V x v is the vorticity 
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FIG. 3: (a) Okubo-Weiss field \(x,y) in the bulk, corresponding to the velocity fields a and b of Figf2] (b) clearly shows the prominent 
saddle (violet) and the surrounding vorticity field (yellow). Note that, in (a) also many saddles (violet) are visible but their intensity is two 
order of magnitude weaker than that in (b). (c) Log-log plot of the distribution of the Okubo-Weiss field. We show separate PDFs' for the 
positive (vorticity) and negative (extensional) values of A, indicated by A+ and A_ respectively (A_ = |A| when A < 0). The PDFs' are nearly 
symmetric and has power law regimes. 
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FIG. 4: Distribution of displacement magnitudes of all the particles, 
collected after different time intervals t (shown in arbitrary units in 
the figure). The inset shows data collapse of the plots shown in after 
rescaling u by t, and P(u) appropriately. This collapse is possibly 
a consequence of the fact that the average displacement (u x ) = 
{v x )t dominates the behavior of u. But the scaling behavior cannot 
be explained by this. The scaling exponent is —1.2 and the plots for 
large t only show significant scaling regimes. 



the strain tensor. Even in viscous flows A turns out to be an 
useful measure and regions with vortices have A > 0, while 
the strain dominated regions have A < 0. Note that a saddle 
correspond to stretching in one direction and compression in 
the orthogonal direction, i.e., it generates a strain dominated 
region. 

Our system is not strictly incompressible, so we computed 
the incompressibility parameter k = ((V.«) 2 )/((V?) 2 ), 
where the denominator is essentially ^ i -(diVj) 2 . For an in- 
compressible fluid the numerator is zero, while our system 
yields a value of k < 0.2 (averaged over many configura- 



tions), which is small enough for an incompressibility approx- 
imation to be valid. Also in terms of total particle numbers, 
the fluctuation is less than 1% (less than 20 in 2000). For com- 
puting k we interpolated the particle velocities onto a square 
grid. In Fig§}a,b we plot the Okubo-Weiss field correspond- 
ing to the velocity fields (a and b) of Figj2] We also compute 
the distribution function of A, shown in Figj3}c. 

Finally we report intriguing power laws in the PDF of the 
particle displacements \Hj\, in the bulk (excluding the bound- 
ary region where (v x (y)) is exponential). Here j is the particle 
index. It turns out that although velocity of the particles, Vj, 
are quite random, the displacements Uj after large time in- 
tervals show characteristic patterns, around the plastic events 
(figure not shown here). The PDF of u = \uj\ is shown in 
Fig. [4] which, at large time intervals t, shows two clear power 
law regimes. Rescaling u with t (and also P(u) appropriately) 
the PDFs' collapse nicely (inset of Fig. 4), although the PDFs' 
for short t do not have any power law regime. 

The different scaling regimes of P(u), namely, oc u and 
oc u^ 1 - 2 , reflect distinct kind of particle motion in the sheared 
polycrystal. P(u) is the fraction of particles undergoing par- 
ticular type of motion and is therefore approximately propor- 
tional to the area fraction occupied by these particles in a typ- 
ical velocity map like Fig. 2b. The displacements are small 
at the core of the large grains where motion is vortical. As- 
suming a slow rotational speed uio, the displacement u, for 
u!ot <C 1, is u ~ uJort, where the radius r is measured with re- 
spect to the center of the grain. Thus P(u)du oc dA = 2irrdr 



and using u ~ ui^rt, we get P{u) 



' (wot) 



Consequently a 



time independent collapse occurs in the P{u)t versus u/t plot 
(inset of Fig.4). Larger displacements (the «~ 1,2 regime) are 
dominated by (v x ) and here approximately (v x (y)) ~ y~ 3 
(log-log plot of (v x (y)) versus y not shown here). Now in 
this case P(u)du oc dA = dy.L and using u ~ y~ 3 t, we 
get tP(u) ~ ^(j) -4 ^ 3 ; again a t independent collapse. Al- 
though our P(u) is restricted to the bulk region, any remnant 
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effect from the boundary region, where (v x (y)) is exponential, 
would contribute a u^ 1 scaling. The observed u~ 12 scaling 
is possibly a mixed effect. 

In summary, we have shown that the plastic flow in 
sheared polycrystals show strong dynamical heterogeneity 
which manifests as three distinct regimes in the displacement 
distribution of the particles. Further, the elementary plas- 
tic events of the flow field can be explained in terms of the 
underlying dislocation dynamics. That bridges two seem- 
ingly disparate descriptions, namely continuum and discrete, 
of sheared solids. 
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